PRECONDITIONERS BASED ON WINDOWED FOURIER FRAMES 
APPLIED TO ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS 
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Abstract. We investigate the application of windowed Fourier frames (WFFs) to the 
numerical solution of partial differential equations, focussing on elliptic equations. The 
action of a partial differential operator (PDO) on a windowed plane wave is close to a 
multiplication, where the multiplication factor is given by the symbol of the PDO evalu- 
ated at the wave number and central position of the windowed plane wave. This can be 
exploited in a preconditioning method for use in iterative inversion. For domains with 
periodic boundary conditions we find that the condition number with the precondition- 
ing becomes bounded and the iteration converges well. For problems with a Dirichlet 
boundary condition, some large and small singular values remain. However the iterative 
inversion still appears to converge well. 



1. Introduction 

Localization in the position-wave number space is an important concept in partial 
differential and other operator equations. The large class of pseudodifferential operators 
acts approximately local in the position-wave number space. Prominent examples where 
this is put to use in numerics are wavelets and multigrid methods, methods that are closely 
related. Multigrid and wavelet approaches are very successful for elliptic problems, and 
have been applied also elsewhere. For wavelet concepts we refer to [7j, the multigrid 
literature is very large, see e.g. the book [19J. 

Wavelets provide a decomposition in scale and position. However, they have little 
resolution in direction, i.e. the direction in the wave-number space. More recently, a 
frame of functions [H] called curvelets has been proposed as a tool for numerical analysis 
of PDEs [3]. Curvelets provide additional localization in direction. The localization 
in direction is better and better for the smaller scales by a so called parabolic scaling: 
A fourfold smaller scale leads to twice better resolution in direction and twice better 
resolution in position. A method of solving a pseudodifferential equation using curvelets, 
including a curvelet based approximation for the operator inverse suitable for use as a 
preconditioner was introducted in [12J. Approximations of pseudodifferential operators 
were discussed in [9]. 

In this paper we analyze windowed Fourier frames, also referred to as Gabor frames. 
They provide a decomposition of the position- wave number space ("phase space") into 
rectangular blocks. Compared to curvelets it is a simpler decomposition, and easier to 
implement. Compared to wavelets it still offers better directional resolution, although 
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at somewhat higher (log-linear) cost. Windowed Fourier frames are also simpler than 
curvelets in the sense that they are generated by translations and modulations of a single 
window function, while for curvelets there is no such set of transformations that exactly 
maps the one to the other. 

Chan et. al [4 J has considered elliptic problems and proposed circulant preconditioners 
to solve the resulting system of equations using iterative techniques, for example, the 
conjugate gradient method. They prove that such preconditioners could be chosen to 
reduce the condition number from 0(n 2 ) to 0(n), where n grid points have been chosen 
to discretize the problem for a second order elliptic problem. Some popular precondi- 
tioning techniques to solve linear systems using iterative methods have been discussed in 
detail in [3], Chap. 1], [11, Chap. 10], [UJ Chap. 7] and references there in. These 
techniques includes use of positive definite matrices, incomplete LU and cholesky factor- 
izations, multilevel, multigrid, wavelet preconditioners and so on. It has been noticed 
that condition numbers are in control with most preconditioners and have slower growth 
compared with the unpreconditioned system. 

We believe that, while multigrid and wavelets are very important as preconditioning 
methods, other possibilities should also be studied. The particular phase space localiza- 
tion of wavelets is well suited for elliptic equations, but not for other types of PDE, for 
example the Helmholtz equation. Windowed Fourier frames, curvelets or maybe other 
transforms are natural candidates to study in more general settings than the elliptic 
problem. 

Our purpose is therefore to introduce a preconditioning method based on windowed 
Fourier frames, and to establish that it performs well for certain discrete PDE's. We start 
with straightforward elliptic problems, and include an example where the different phase 
space localization properties provide an advantage for windowed Fourier preconditioning. 
For this study we consider a symmetric second order elliptic BVP with a finite and a 
periodic domains. We first discretize the PDE using a standard finite difference scheme. 
We use a preconditioner based on windowed Fourier frames (WFFs) and the symbol of 
the operator while solving the discrete PDE using iterative linear solvers, to speed up the 
convergence. 

The article is organized in the following way. We start by introducing the precondi- 
tioners in Section [2| We study boundedness and invertibility of a symmetrically precon- 
ditioned operator in Section [3j Some numerical test results are presented in Section [4j 
We finish this study in Section [5j with some concluding remarks. 

2. Preconditioners based on the symbol of the operator and a 

windowed Fourier frame 

In this section, we focus on introducing and defining the preconditioner based on the 
symbol of the partial differential operator and a windowed Fourier frame. As a model 
problem we consider the boundary value problem 

(1) Cu = f in n, 

with Dirichlet and periodic boundary conditions on dfl where 

Cu = — V • (a(x)Vu) + b(x)u(x), 
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for given real coefficients a(x), and b(x). Here we consider Q C M d , an open bounded 
domain, / : Q — > R is a given function and u : Q U <9fi — )• R is an unknown func- 
tion. A detailed description of this type of problems can be found in [21 QUI ITB] . and 
many references therein. We investigate both periodic and non-periodic boundary value 
problems. 

We make the following assumptions on the coefficients, to ensure that £ is boundedly 
invertible. We assume there is some constant C such that a(x) > C. In the case of Dirich- 
let boundary conditions we consider b > 0, in the case of periodic boundary conditions 
we assume that b > Co > for some constant Co. 

Windowed Fourier frames. We first give a short introduction of windowed Fourier 
frames (WFF). Let H be a Hilbert space. A sequence {^ n }ner is a frame [TH Section 
5.1.1, Definition 5.1] of H if there exist two constants A > 0, B > such that for any 
/ 6 H, 

(2) A\\f\\ 2 <J2\(fM\ 2 <B\\f\\ 2 . 

ner 

The index set T might be finite or infinite and one can define a frame operator Fx so that 

F 1 f[n} = {fM, for all n G T. 

If the condition (|2| is satisfied then Fx is called a frame operator. When A = B the frame 
is said to be tight [HJ page 155], [6]. 

A window function is simply a function in C oo (R <i ), i.e. it is smooth, and zero outside 
some chosen finite domain. Let us focus on d — 1, and consider a real symmetric (g(t) = 
g(—t), for all t E 1), non-negative window function. We can translate g by v G K and 
modulate g by frequency £ G K as g v ^{t) = e l &g(t — v), which are known as windowed 
Fourier atoms or Gabor atoms. Here = 1 for any v G R, £ G R since ||g|| = 1. 

If the functions g v ,^{x) satisfy the frame condition (|2]), then they are called windowed 
Fourier frames (WFF) [TJ1 Section 5.4] and for any / G L 2 (IR), the operator F defined 
by 

(3) Ff(v, = (/, g Vt z) = [ f(t)g(t - v)e~^dt 1 

is called the windowed Fourier frame operator. In Figure [T] we present a sample windowed 
Fourier frame. 




Figure 1. Example of a function from a windowed Fourier frame in one 
dimension with window function w(x) = sin e -c/x e + ^-c/(i-x) J , c — 1-5. 
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A version of the windowed Fourier frame transformation with discrete index set T, is 
obtained by restricting (v, £) to a rectangular grid with interval size v o and £o i n time 
and frequency respectively and define [HJ page 182] 

9n,k = 9(t ~ nv )e^ k \ 

where = k£o, which will be needed in Section [4] while implementing the precondi- 
tioners for elliptic PDE's. It is well understood 0, [H] that the windowed Fourier fam- 
ily {g n ,k} (n,k)ez 2 is a frame only if > 1 and the frame bounds A and B satisfies 
A < < B. Let g be a window function with support [—^, £]■ Then {g n ,k}r n k)^z 2 * s 
a tight frame [H] with a frame bound equal to A, if 

(4) — V \g(t - nv )\ 2 = A > 

« n=-oo 

for all t e R. 

In the discrete setting, we replace the Fourier basis {e lk ^ ot }kez of L 2 [— 7r/£ , 7r/£ ] by 

i27rkn 

the discrete Fourier basis {e ^ }o<k<N t of C ; to construct a discrete windowed Fourier 
frame. Let us consider g[n] be a N periodic discrete window function with a support and 
restricted to [-N/2,N/2] that is included in [-JVj/2, JVj/2 - 1]. Then following [H] one 
can see that if mod (N, M) = and 

M x 

Ni \g[n-mN}\ 2 = A > V < n < N 

m=0 

then {g m ,k[ n } = g\ n ~ m M}e t2nkn ^ Nl } < : k<N l ,o<rn<K is a tight frame in with frame 
bound equal to A where K = jj. For a fixed window position indexed by m, the discrete 
windowed Fourier frame coefficients are 

N-l 

i27vkn 

f[n]g[n - mM]e »t , 

n=0 

for all < k < N t . 

We set the length of the window equal to 2v , so that two successive windows have 
overlapping support and we choose tight windowed Fourier frames. So we choose the 
parameters such that |^ = 2v , and g has support in [— v ,Vq\. 

Several choices for the window function are discussed in Appendix |A} Our initial crite- 
rion for a good window function is the decay of its Fourier transform. We conclude that 
both the windows h^(x) with c = 1.5, d = 0.9 and the window formed by h±(x) are very 
well behaved. We experimented with both the windows to define the preconditioners and 
recorded the number of iterations needed for convergence while solving the precondition 
linear system. The numbers for these two choices of the window function are approxi- 
mately equal, preconditioners formed by using the stretched window h*>{x) with c = 1.5, 
d = 0.9 take a few iterations less than that of h^x), but the effect is small. We decide to 
use the window h 5 (x) for this study since we designed it (and it performs a little better). 

In order to apply windowed Fourier frames to the BVP ([!]) we need to define and 
organise window functions for a multi-dimensional and bounded domain Q C M d , say 
Cl = [a, b] d , where the domain can be both periodic and non-periodic. Multidimensional 
window functions are defined straightforwardly using a tensor product approach. For 



periodic domains we assume that the period is a discrete multiple of t>o, say Kv . Then 
the set of coefficients becomes periodic with period K . In Section [3] we will discuss this 
further. For non-periodic problems, we define K + 1 windows in each of the coordinate 
directions (instead of K windows) with support ^ (in each of the coordinate directions), 

l c 

K ' " 1 K ■ 



considering an extended domain [a — b + D f2, by 



9A X ) 



g(x - jv ) iGfi, 

o x i n, 



where j G {1, 2, • • • , (K + l) d }. Then we have 2 d (K + 1), (d > 1) windows that cross the 
boundary of f2. 

Preconditioners. Here we return to our main discussion. First we give a short mo- 
tivation why a PDO can be approximated by WFF's and the symbol of the PDO. For 
simplicity we start with BVP (HI) in one dimension. Setting 

u = 9j,k(x) = e iikX g{x - jv ), xeQcR, 

j £ J (the set J is defined in Section [3]), k G Z, we get 

d^u d 
a{x) dx^ = a{x) Yx 



'G,< ' //(•'' - jv ) + e ^ kX ^9(x - jv ) 
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= a(x) '//(./• - jv ) + (1 + ^k)e^ x ^g(x - jv ) + ^ kX ^g(x - jv 0/ 

= A 1 + B + C, 

where A\ is the leading part (when is very large) of the elliptic operator acting on 
([!]). When ^ is very large, £f is the dominating term in a(x)J^ and so the differential 
operator a(x)-^ can be approximated by a multiplication operator a(x)£%. Now let us 
represent a function m(x) G H using windowed Fourier frames 

where 

Cjjt = (u(x),g jjk (x)} and g jjk {x) = e^ kX g(x - jv ). 
Then we can approximate 

a ( x )f~^ ~ ^2 a ( x )£l C hk9j,k( x ) = "£2 a ( x )£k( u i9j,k}9j,k- 

Here we observe that the PDO (a(x)|^) can be approximately reconstructed by the tight 
windowed Fourier frames and its duals multiplied by the symbol of the PDO if a varies 
little over the support of g^. This result motivates us to define preconditioners based on 
symbols of the PDO's and WFF's. 

It is to note that we consider Xj, j G J, as midpoints of the subdomains of Q. To 
bound the condition number and to speed up the convergence of iterative solvers we are 
interested in defining preconditioners for a discrete equivalent of ([I]) based on the symbol 
and WFFs. To this end, we define an invertible matrix M 

M ikj,k = 5 iA4 1 + a ( x ^- 
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We solve an equivalent system (of Au = f, in which A is a symmetric discrete equiv- 
alent of the operator acting on ([I]) using a finite difference scheme). Several forms of 
preconditioning can be distinguished: 

Symmetric Preconditioning: PAPu = Pf, coupled with Pu = u, where P = 
F*M •/•. 

Left Preconditioning: PAu = Pf, where P = F*M~ 1 F. 

Right Preconditioning: APu = f, coupled with u = Pu where P = F*M~ l F. 

In each case F* is the conjugate transpose of the frame operator F. The relation between 
the left and right preconditioners is established in the following result. 

Theorem 2.1. The singular values of the left preconditioned matrix PA and the right 
preconditioned matrix AP are equal where P = F*M~ 1 F. 

Proof. We have 

(PA)* =A*P*= AP, 

since A is symmetric and 

P* = (F*M~ 1 F)* = (M~ l F)*F = F*M F. 
proves the claim since singular values of PA and (PA)* are the same. □ 

3. BOUNDEDNESS AND INVERTIBILITY OF THE SYMMETRICALLY PRECONDITIONED 

OPERATOR PCP 

A key property of a preconditioner, is that the condition number of the preconditioned 
operator is bounded independent of the discretisation. We will show that the symmet- 
rically preconditioned operator PCP, where P = F* M~ X I 2 F and F stands for tight 
windowed Fourier frame operator, considered in the continuous case, defines a continu- 
ous map, with continuous inverse on L 2 (Q). Under a suitable discretisation this can be 
used to derive the first mentioned property, although we will not do so (instead we study 
some discrete systems numerically in the next section). We prove this in the case of the 
domain 

Q, = [0, 2vr) n = T n 

i.e. [0, 2ir} n with the periodic boundary conditions. In the case of a domain with a 
boundary we will see in the section on numerics that there are a few singular values that 
grow if the discretisation parameter becomes small, so that the result appears to be false. 

Theorem 3.1. The self adjoint operator PCP : -L^(^) — > ^(Vt) is boundedly invertible, 
that is, there exists ci > c\ > such that 

c i|M| 2 < (PCPu,u) < C2||w|| 2 , for all u E H. 

By assumption C : H 1 — > H^ 1 is boundedly invertible. Therefore, to establish The- 
orem 3.1, it is sufficient to show that P is boundedly invertible from L>2(fl) to H 1 ^), 
and that P is boundedly invertible from if~ 1 (fi) to L 2 (fi). In fact it sufficient to show 
that P is boundedly invertible from L-2(VL) to H l (VL), as this implies that P* is boundedly 
invertible from H' 1 ^) to L 2 (0) and P = P*. 

One may think of T n as the hypercube [0,27r[ n C M n or [— 7r,7r[ n . Functions on T n 

may be thought as those functions on IR n that are 2n periodic in each of the coordinate 
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directions. Let p G R. The sobolev space H P (Q) is the space of all functions if) G L 2 (Q) 
that satisfy 



(5) ing = E( 1 + h 2 ) p 



dm < OO 



for the Fourier coefficients a m of ip. The space H P (Q) is a Hilbert space with the scalar 
product defined by 

(6) (0,W= ^ (1 + |m| 2 ) p a m 6 m 

for <fi, ijj G H P (Q) with Fourier coefficients a m and 6 m , respectively 

When we consider one window only, the windowed Fourier frame operator F becomes 
the Fourier transform operator. So we first prove the invertibility of the operator PCP 
considering F as the Fourier transform operator. From here we denote Mi = M~% . We 
first show that Pu G H l {&) for u G L 2 (0). Here 

(Pu,Pu) = (F^MiFu.F^MiFu) = (MiFu,MiFu 



where Fu G £ 2 (Z) is the sequence of Fourier coefficients of u G L 2 (VL). Now using ^ we 
have 



(Pu,Pu) 1 = (1 + |m| 2 )(Mi) 2 |(Fw 



I I 2 

'ml j 



where (Fu) m G ^ 2 (Z) are the Fourier coefficients of u G L 2 (f2). Now for any m G Z n 
there exists < A < B such that A < (1 + |m| 2 )(Mi) 2 < B. So one gets 



2 ' 

12 ^ II D„.ll ^ DIL.II2 



< \\Pu\\ x < B\\u\ 
and 

Pu = F^MiFue H\n), ifuGL 2 (0). 
Thus there exists < c\ < c 2 such that 

c i|M| 2 < (PCPu,u) < c 2 ||u|| 2 , 

and that confirms the boundedly invertibility of PCP : L 2 (fi) — > L 2 (f2), since C is 
boundedly invertible and P* = P. 

Next we investigate the boundedly invertibility of PCP : L 2 (f2) — >■ L 2 (fi) considering 
F as a windowed Fourier frame operator. To define window functions we consider K 
uniform subintervals of size ^ in each of the coordinate directions so that we can divide 
Q (= T n ) into K n subdomains. Then we define K n overlapping subdomains with length 
^ on each of the coordinate directions and denote them by Dj C f2 for all j G J where 
J = {1,2, ••• ,K} with K = K n . We consider K windows on each of the coordinate 
directions. We define window functions gj(x) on Q for all j G J satisfying 9]( x ) = 1> 
where x G Q. The support of each gj(x) is contained in Dj. By Xj we denote the midpoint 
of the block Dj. Of course for the periodic setting one needs to consider appropriate 
boundary domains, in one dimension for example, D\ = 2n [0, U 2ix [^^, l). We will 



also assume that the support of gj stays away some small distance from the boundaries. 
The windowing operator is denoted by W, defined by 

(7) W : L 2 (Vl) — ► L 2 (Di) x L 2 (D 2 ) x • ■ ■ x L 2 (D R ) 

(8) Wu = {g x u,g 2 u, ■ ■ ■ } g k u) . 

(Wu can also be viewed as an element of the larger space L 2 (Q,C K ).) The adjoint 
windowing operator W* then equals 

W*(ux,u 2 ,u 3 , ■■■ ,Ux) = ^gjUj. 

The Fourier transform of the windowed Fourier transform is done on each domain Dj. 
We will write Tjy for this Fourier transform on functions restricted to Dj. We redefine 
the operator P as 

We first study P as an operator H l l 2 — > H^ 1 ^ 2 , starting with a result about W. 
Lemma 3.2. Let W be defined by ([SU. Then there exists < a < (5 such that 

(9) a||Wu|| JT i/a(n,c*) ^ \\ u \\h^(si,c) < PWu\\ H y 2{flfiit) . 
Proof. We consider separately the two inequalities 

(10) a\\Wu\\ Hl/2{ ^ CIi) < ||«||Hi/a ( n,o 
and 

(11) IMU*(n,o < P\\Wu\\ 



Equation ( 10 ) is easy, it follows directly from the continuity of the multiplication by gj 
on H 1/2 . Equation rtllk is more difficult. We solve it as follows. Define Er s \ to be the 
pseudodifferential on £1 (= T n , the torus) with symbol (1 + ||^|| 2 ) s ^ 2 - This defines a self 
adjoint operator. The H s norm of a function u is equivalent to ||£J( s )ii||x, 2 . In this proof 
we in fact use this formula for the norm. 
We start with the basic estimate 

(12) \\Wu\\ 2 Hl/2 > l\\Wu\\l 2 + ^\\E (1/2) Wu\\ 2 

where C > 2 is a constant to be chosen later. For ||Wu|||, we find the following 



iwilia = (u,^2g 2 u) = ^2(gjU,gju) = \\Wu\\l 2 . 



j 



For ||i?( 1 /2)W^M|| 2 we use the following 

\\ E (i/2)Wu\\l 2 = ^2(E {1/2)gj u, E il/2)gj u) 



(E {1/2) u,E {1/2) ^g 2 u) + ^(u, [E {1) ,gj]gju}, 
j j 



\ U \\ H l/2 



3 
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where we used that i?( 1 / 2 )-E'(i/2) = ^(i) and [E^,gj] = E^gj — QjE^y The operator 
gj] is a pseudodifferential operator of order (follows from [TJ). Therefore we find 
that there is a constant D such that 



<d\H1 2 . 



We now use what we just described in (12) and we choose C = min(2D, 2). This yields 

1 

in/./ , , 1 1 - / ) 1 1 , , i 

\L 2 j 



\Wu\\ 2 Hl/2 



> 



1 



> 



2 
1 

C 



\u 



\* _i 

\La T C 



\ E {l/2)U\\l 2 



\E(1/2)U 

1 



2 

!l 2 



C" 



2 



We therefore have proved (11). 



□ 



Using the lemma we can prove the following 

Proposition 3.3. P is boundedly invertible if _1 / 2 (f2) — > H 1 / 2 ^). 

Proof. The operator W maps from i? 1//2 (f2) to the Sobolev space H^ 2 (Di) x H 1 ^ 2 (D 2 ) x 
• • • x if 1 / 2 (D^). The proposition follows from the lemma we just proved and the Lax- 
Milgram theorem, using the Fourier coefficient description of this Sobolev spaces. □ 

Next we study P as an operator L 2 (Q) — > iJ _1 (fi), by first proving that P is a pseudo- 
differential operator. 

Proposition 3.4. If the window functions gj(x) G C£°(Dj), then the operator P defined 
by 

is a periodic elliptic PsDO of order —1 where Sj = T^M^^Td, is a periodic elliptic 
PsDO of order —1. 

Before we go to prove the proposition, we will need some definitions and results about 
periodic pseudodifferential operators [TBI 123] • Let V(T n ) be the vector space C°°(T n ) 
endowed with the usual test function topology. Then any continuous linear operator 
A : T>(T n ) — > T>(T n ) can be represented as 

(13) (Au)(x)=J2°A(x,0u(0e lx< , 
where 

(14) a A (x,0 = e~ lx< Ae ix< 

is called the symbol of the operator A. Let m e K and < 5 < p < 1. In our case 5 = 
and p = 1. An operator defined by (13) is called a periodic pseudo-differential operator 



(PsDO) of order a if the unique function a A G C°°(T n x Z n ) defined by (14) satisfies 
(15) \A«d?a A (x,0\ < C^ m (0 m - plal+m 

for every x G T n , for every a, (3 G N n and (f ) := (1 + |£| 2 ) 1/2 . Here by S™ s (T n x Z n ) we 

denote the space of functions o A G C°°(T n x Z n ) that satisfies pHb. If a A G >S™ 5 (T n x Z n ), 

one may denote A G Op S™ 5 (T n x Z n ). An operator A G Op S™^T n x Z n ) is called elliptic 
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if, in addition, the symbol a A G S™ 5 (T n x Z n ) of A satisfies 01 > C \€\ m where 

x G R n , C > and |f | > £ > 0. 



Proof of proposition 3.4 First we show that the operator Sj = J r ^M 1 / 2 J r D . is a periodic 
elliptic PsDO with period ^ on each of the coordinate directions. Now for any u G D(T n ) 
we have SjU = J 7 ^ 1 (M^ (^DjU^))) , with symbol a 3 -(£) = Mi(£) = J — Consider 

ffe G Z n with (ffc)fc = 1 and (t>fc)j = if i ^ /c. Now using difference calculus 

Aj fc Mi(0 = Mi (e + v k ) - Mi(0 = (1 + a(x,-)« + «*)T 1/a - (1 + a(^-)£T 1/2 , 
gives 

i A i M i(Oi<^(i + iei)- 2 , 

and similarly for higher order differences, and thus Mi(£) G S"~](T n x Z n ) [151 [20] . That 
is to say that Sj is a periodic elliptic PsDO of order m — — 1 with period -5 on each 
of the coordinate directions (using the binomial theorem and the Leibnitz formula for 
differences [H] one can deduce the similar relation for any a). The previous statement 
implies that the distribution kernel of Sj, let us denote it here by K(x, y), has singularities 
at x = y + j^rk, where k G Z n . Therefore, when it is restricted to x G supp y G 
supp(^j), then the singularities are contained in the set x = y. 

Then we show that Pj = gj(x)Sj(x, £)gj(x) is a PsDO for all j G J. Here we consider 
gj(x) G C^(Dj), j G J, and thus gj is a PsDO of order 0. Then it follows from [T31 |5Tj 
that the composition Sjgj is a periodic elliptic PsDO of order —1, and so is the operator 
(since Sj is a PsDO of order —1) Pj with principal symbol crj(x,£) = gj(x)cij(£)gj(x). 
Since gj(x) G C™(Dj), it can be viewed ) G C£°(f2). Thus for all j G J the symbols 

<Tj(x,£) can be extended to cr, G S*"]^ x Z n ). 

We have P = Ylj^j 9jSj9j- Here gjSjgj, j G J are periodic elliptic PsDOs of order —1. 
Let us denote the principal symbol of Sj by a,. So it follows that P is a periodic elliptic 
PsDO of order —1 since P is a sum of PsDOs of order —1. □ 

Since P is an elliptic PsDO, it follows from [T31I2T] and references there in that there is 
a parametrix Q such that QP = I+R with R a smoothing operator. For our purposes it is 
sufficient that R is a PsDO of order —1 (see PsDO literature, parametrix, for example [TJ 
page 18]). Let us now discuss the kernel of the operator P. The kernel of P is contained 
in the kernel of QP. The operator R is compact. We work with operators on L 2 , so it 
is an operator from L 2 to H 1 , and therefore compact as an operator on L 2 . But since 
this is PsDO theory, we can also work on operators on H s , and then R is continuous 
H s — > H s+1 and compact as an operator on H s . From the theory of compact operators, 
it follows that, for any C > 0, there is at most a finite dimensional set of vectors in L 2 , 
such that ||i?u||i 2 > C||w|| (e.g. set C = 1/2.) Then it shows also that the kernel of PQ 
is finite dimensional and hence that the kernel of P is finite dimensional. Similarly the 
cokernel of P, which is defined as the kernel of P is finite dimensional. We claim that 
the elements of ker(QP) are in C°°. The elements of ker(QP) satisfy u = —Ru. Let us 
say u G L 2 . Then Ru G H 1 because R is a PsDO of order —1. But u = —Ru, so in 
fact u G H 1 . But then Ru G H 2 , because R is continuous H 1 — > H 2 . And so forth. 



It follows that u G C°°. Since we already proved in Proposition ^3 that P is invertible 
H~ 1 / 2 — > H 1 / 2 , it follows that the kernel and cokernel of P are the zero sets and that P 
is invertible H s H s+1 , set 
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Next we study the inverse P 1 and show that it is a PsDO. As P is an elliptic PsDO, 
there exists Q G Op(5 1 ) [131 EI] such that QP = I + R and PQ = I + R, where 
R G Op(S-°°) and R G Op(5 ,_0 °). Let us consider V : P 1/2 H- 1 ' 2 be the inverse of 
P : H- 1 ' 2 ->■ PT 1 / 2 . Then one gets 

(16) qpv = (/ + p)v 

and 

(17) QPV = Q. 



Combining (16) and (17) we have (J + R)V = Q and so 

(18) Q - V = RV. 
Also one can write 

(19) VPQ = Q 
and 

(20) VPQ = V(I + R). 



Combining (19), (20) and (18) one gets 

Q — V — VR = (Q + (V - Q))R = QR- RVR e OP(S~°°), 
and so V G OP(5 1 ). 

From here we may conclude that Theorem 3 . 1 1 is proved. However, there is no estimate 
of Ci, so we add a slightly more direct estimate for the constant c\. 

Proposition 3.5. There are constants < A < B such that 

(21) A||Pu||_g-i < ||m||l 2 < -B||Pw||#i, f° r an U u G L 2 (Q). 

Proof. The second inequality follows from the fact that P is a pseudodifferential operator 
of order 1. For the first inequality consider the family of pseudodifferential operators 

/ 2 \ s/2 

E\ <s , with symbol E\, s (x, £) = ( 1 + |a ) • Clearly Pm satisfies 

> \\Ex tl / 2 Pu\\ H i/2 = \\PE x ,i/2U+ [E x ,i/2,P]u\\ Hl/2 

\ — — J II || || || 

> \\PE XA/2 u\\ Hl/2 - ||[P Ai i /2 ,P]m|| h1/2 . 
The first term can be estimated by 

,1/2^11^1/2 > V^||-f > ||ii"- 1 /2_»i2-l/2 |M|l 2 - 

The second term can be written as 

( 23 ) X}[-Ea,i/2, = ^[^,1/2,^]^^ + gjSjlE^, gA 

3 3 

The commutator [E Xt i/ 2 ,gj] is a pseudodifferential operator whose symbol can estimated 
in terms of A, it is given by an asymptotic sum of the form 

oo 
3=0 

where the Bj are symbols of order —1/2 — j. It follows that we have the estimate 

(24) \\[E x , 1/2 , gj }\\ H ^ Hs <C\- 1 . 
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Equations (23) and (24) show that the second term on the r.h.s. of (22) can be estimated 
by CX^ 1 \\u\\ L 2. By choosing A large enough, it follows that the first inequality of (21) 
holds. □ 



4. Numerical experiments 

Now the question arises how efficient it is to solve constant and variable coefficient 
PDEs using the proposed preconditioner? Does the preconditioned system converge faster 
than the unpreconditioned system while using some iterative solvers for linear system of 
equations, for example, the conjugate gradient method? To investigate it, we experiment 
our technique focusing elliptic BVPs and compare it with the conjugate gradient method 
within the framework of Matlab. 

4.1. One dimensional examples. Let us start with a one dimensional case (d = 1). 
Then the BVP ([I]) takes the form 

(25) - 4- K*)^r^ ] + b(x)u(x) = f{x), V xeQcR 



dot; V dec 



with 



u = on dQ 



where a(x), and b(x) are known functions of x with a(x) ^ 0. We approximate the 
problem Cu = f using a finite difference scheme. Let us define the grids as Xj = jh, 
j = 0, 1, 2, ■ ■ ■ ,N with grid spacing h, and N > 2. We define the approximations to the 



solution u by Uj defined on the grid points Xj for all j = 0, 1, 2, • • • , N. The BVP (25) 
can be approximated by 

(26) - — (a i+ xu i+1 - Ui(a i+ i + a^i) + a^xm-^j + hui = f, 

for all i = 1, 2, - ■ • ,N where a { = a(xi), u { = u(xi), b { = b(xi), f = f(xi), a i± i = ct ' + °' ±1 



and using boundary condition u = = 0. We write (26) in the matrix form as Au = f. 
In Figure |2j we plot spectral radius and condition numbers of A for several choices of 
system size to demonstrate the polynomial growth of condition numbers. It is noticed 
from Figure [2] (also several articles and books on finite difference schemes have a clear 
discussion about the eigenvalues of A, e.g., [21 [15], [TF] and references therein) that the 
condition number of a matrix A grows like 0(1/ h 2 ) (when a(x) = 1), and as a result the 
convergence of any iterative methods become slower. 

We consider three examples to address the questions arise in our previous discussion. 
In the first example, we consider a(x) = 1, b(x) = and observe singular values and 
condition numbers of the symmetric and left preconditioned matrices. We also compare 
the number of iterations taken by the preconditioned solvers to converge. Then we 
repeat our experiments when a(x) varies with x in the second example. Here we compare 
condition numbers by varying number of windows (to cover the domain) to show the 
efficiency of windowing. We also show the advantage of using the preconditioned CG 
solver by comparing the number of iterations with the CG solver. In the third example 
the coefficient varies discontinuously. 
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Figure 2. We show spectral radius, and condition number of the matrix 



operator acting on ( 26 ) for various choices of system sizes with a = 1 and 
6 = 0. 



Example 4.1. Consider the BVP 

d 2 



dx' 



u{x) = f(x), V < x < 1, 



with boundary conditions 

(1) u{0) = and u{l) = 

(2) periodic boundary condition u(0) 



u(l) 



The discrete operator A can be found in (26), to enforce the periodicity, when needed, 
we define A(N,1) = A(l, N) = —j^- In Figure [3] we compare singular values of PAP 
and PA with the singular values of A (with periodic and non-periodic settings). From 
these computations, we notice that most of the singular values are close to one, whereas 
the singular values / eigenvalues of A are [13 page 561] 



A, 



sin 



JIT 



V J = 1,2, 



N. 



h 2 \2(N + 1) 

In the periodic case there is a singular value zero due to the constant solutions. 

In Figure [4] and Figure |5j we show the condition numbers of the symmetric and the 
left preconditioned operators considering A as a periodic and a non-periodic operator re- 
spectively. We notice that the condition numbers of PAP and PA appear to be bounded 
when A is periodic. The condition number of PAP grows as O(N) whereas condition 
number of PA grows similar to the unpreconditioned operator A when A is non-periodic. 



We also notice from Figure 



that 



a 3 



grows much slower or hardly at all for both the 



Fourier frame preconditioned operators PAP and PA (non-periodic operators), whereas 
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Figure 3. Here we plot singular values of the preconditioned operators 
(non-periodic operator (left figure) and periodic operator (right figure)) 
where K = 4 and iVj = 2 6 . From here A appearing on the figures represents 
the discrete operator A. 



grows as of ^ for the unpreconditioned operator A (A3 is the third highest singular 



a 3 

value and \n~2 is the third lowest singular value). Thus analyzing Figure |3j Figure |4j and 
Figure [5] we expect that the left preconditioner and the symmetric preconditioner would 
perform well since most of the singular values are clustered to a constant for the both 
preconditioned system (though the condition number of the left preconditioner grows as 
of A when we consider a non periodic boundary condition). The condition number of 
PA and AP are the same, so we present experimental results of the left preconditioned 
system only. 

Then we focus on solving the problem using the conjugate gradient solvers. We present 
the number of iterations taken to converge for three different choices of solvers considering 
f(x) = exp(27r(a; — 0.5)) in Figure[6] From this experiment we notice that both the SPCG 
and the LBICG take a very few iterations to converge compared to the CG (for both the 
periodic and the non-periodic boundary value problems). 

This is to note that we have considered windowed sine frames (WSF) as well for the 
same computations. It produces the similar results as of WFF, so one can also use WSF 
for the computations as well. 

We have experimented with the preconditioner considering a simple constant coefficient 



BVP in Example 4.1 This preconditioner is actually designed for variable coefficient 
problems. In the following example we aim to demonstrate the advantage of windowing 
for such problems. 



Example 4.2. Consider the periodic BVP 

d / du(a 
dx \ dx 



' a{x)^y^- ) + b(x)u{x) = f(x), V < x < 1, 
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Figure 4. Condition numbers of the operators for various choices of sys- 
tem sizes (periodic case). 
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Figure 5. Condition numbers of the windowed Fourier frame precondi- 
tioned operators (non-periodic case) for various choices of system size. 



where a(x) = 10 — 9.5 cos(27rx), b(x) = 1 and 

, . J exp(x) when < x < 0.25 
= \ exp(-x) if 0.25 < x < 1. 

We use the scheme (26) and consider tol = 10 -10 . In Figure [7J we compare condition 
number of the discrete symmetric preconditioned variable coefficient operator in < x < 
1 by varying the number of windows to cover the domain. From this experiment we notice 
that the more subdivisions of the domain (using windows) one considers, the smaller the 
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FIGURE 6. Number of iterations taken by various solvers to converge (we 
consider 4 windows here) considering boundary conditions (1) u(0) = 
u(l) = (left figure), (2) u(0) = u(l) (right figure). 

condition number becomes. Then we compare the number of iterations taken by the 
CG method and the SPCG method (considering 8 windows to cover the periodic [0, 1] 
domain) and notice the good behavior of the preconditioned system. 
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Figure 7. We show the condition numbers for various choices of use of 
number of window functions. This left figure shows the efficiency of win- 
dowing by reducing the condition number. The right figure shows the 
efficiency of windowing by reducing the number of iterations significantly. 



Example 4.3. Consider the BVP 

d ( , s du(x 



dx V dx 



+ b(x)u(x) = f(x), V < x < 1, 
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with boundary conditions u(0) = and u(l) = 0, where b(x) = 1 for all < x < 1, f(x) 
as a random function and 

exp(x) when < x < 0.5 
exp(— x) if 0.5 < x < 1. 

We use the scheme (26 ) and consider tol = 10 -10 . In Figure M we notice that the condition 



a(x) 



number keeps growing but more slowly for large N, while the number of iterations appears 
to becomes constant or close to constant. 
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Figure 8. The left figure shows condition numbers for different choices of 
windows, and the right figure shows the number of iterations taken to con- 
verge by the conjugate gradient method and the symmetric preconditioned 
conjugate gradient method. 



4.2. Two dimensional examples. Here we apply our technique in two dimensions. We 
consider one example where the operator is very strong in one direction to demonstrate 
the advantage of using the exact symbol. For the other example we consider an operator 
that contains strongly varying coefficient a(x,y) acting on it to demonstrate further the 
advantage of windowing. 



Example 4.4. Consider 
Cu(x,y) 



-10 



d 2 u(x,y) 1 d 2 u(x,y) 



f(x,y), 



dx 2 10 dy 2 

and Q C IR 2 with u\qq = 0. Now let Q = (0, l) 2 be the unit square and similar 
to one dimensional case, we define = Va; = h = 1/N > and Xij = (ih,jh), 
i, j = 0, 1, 2, • • • , N. The results have been compared in Figure [9] Here we notice that 
the windowed Fourier frame preconditioners designed with the exact symbol performs 
significantly better. 

Example 4.5. Consider 



Cu(x,y) 



d_ 

dx 



a\x) 



_ du(x, y) 
dx 



d_ 
dy 



Ky) 



du(x, y) 
dy 



f{x,y), 
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Figure 9. In the left figure, we compare number of iterations consider- 
ing the exact symbol S(£, rj) — 1 + 10£ 2 + ^rf and an isotopic symbol 
S(£,r]) = 1 + £ 2 + rf. Here we use N l>x = N hy = 2 2:7 (thus N = 2 5:10 is 
the number of points on each direction). In the right figure we compare 
relative residuals for each iteration number considering Ni jX = N^ y = 2 4 . 
For this computation we consider K x = K y = 4 windows in each directions, 
f(x,y) = e-t* 4 *) and tol = lO" 10 . 

in a periodic domain Q = [0, l] 2 , where a(x) = 10 — 9.5 cos(27rx) and b(y) = 1. For this 
computations we consider symmetric preconditioned solvers only. In Figure [TUJ we com- 
pare the Fourier preconditioned conjugate gradient method (one Fourier transform only, 
considering no windowing) and the windowed Fourier frame preconditioned conjugate 
gradient method (4 windows in each direction). Here we observe that the preconditioned 
solver based on the windowed Fourier frame performs better than the Fourier precondi- 
tioned solver. 

5. Conclusions 

We study windowed Fourier frames (WFF) focusing on elliptic boundary value prob- 
lems and present new preconditioners based on the symbol of the operator and WFFs. 
From this study we conclude that the preconditioners based on WFFs work well for el- 
liptic problems, most of the singular values clustered to a constant approximately for 
both preconditioned operators (SPCG and LBICG). For periodic domains, the condition 
number is bounded. However, when a non-periodic domain is considered, there are a few 
very large and a few very small singular values, which cause the condition number of 
the symmetrically preconditioned system to grow O(N) whereas the condition number of 
left preconditioned system grows close to 0(N 2 ). This affects the convergence of the CG 
method relatively little, as it concerns only a few singular values, but a better definition 
of such a preconditioner on a bounded domain is clearly a question for further research. 

As expected both the SPCG and the LBICG take very few iterations to converge com- 
pared to the unpreconditioned CG method. For multidimensional PDEs, if the problem is 
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Figure 10. We compare number of iterations considering symmetric pre- 
conditioner with 4 windows and without windows. We consider the symbol 
V ) = 1 + a(x)e + b(y) V 2 . Here we use K x = K y = 4, N hx = N l>y = 2 3:7 
(thus N = 2 6:10 is the number of points on each directions). For this com- 
putation we consider f(x,y) = e~^ x+v ^ and tol = 1CT 10 . 



strongly dominant in one direction then the use preconditioners based on the exact sym- 
bols is recommended. We notice that the use of reasonably many windows (while defining 
the preconditioner) has an advantage for PDEs with strongly varying coefficients. 



Appendix A. Window function construction 

The window function plays an important role in the preconditioner, therefore we discuss 
its construction. A priori, the objective is to select a window function that satisfies the 
properties described in section [2] and decays rapidly in the frequency domain. It is well 
known that smooth functions have fast decay for large frequencies [14] . So we want to 
form window functions that have at least some vanishing derivatives at both ends of the 
support. In [14] there are some general rules for constructing window functions, as well as 
some specific examples, we make use of this and also compare some of the example with 
a construction of our own. We present a schematic window function and two translations 
in Figure [11] 

To design g(x), we consider a monotone increasing function h(x) (C k or C°°) such that 



(27) h(x) 
and satisfies 



if x < -1, 

1 if x > 0, 



(28) h 2 (-l/2 + y) + h 2 (-1/2 -y) = 1, -1 < x < 0. 
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Figure 11. The figure shows window function g(x) has support [—1,1], 
the left and the right half graphs are the translations of the function. 

Then to form a function g with support on [— v ,v ] we take 

J h(x/v ) for x < 
^ \ h(—x/vo) for x > 0. 



then (28) ensures that Q holds. Of course a variation is possible. For example we can 
take a parameter 1/2 < a < 1, and, let g(x) = for \x\ > a, g(x) = 1 for |x| < 1 — a, and 
g described by a dilated version of h for 1 — a < \x\ < a. 

In Figure [12J we show some examples of monotone increasing functions. First we have 
a class of functions, with different c, given by 

{0 if x < 0, 

sin (% c e ~t a ) if 0<x<l, 
1 if x > 1, 

(we need to translate this backward by 1 to have a function /i satisfying the description 

_ i 

above). This is constructed using the standard profile h a (x) = _ , e x r ~ , that goes 
smoothly from to 1 on [0, 1] and satisfies h a (x) + h a (l — x) = 1. Taking the sin(|/i a (a;)) 



instead of h a transforms the property h a (x) + h a (l — x) = 1 into the similar property (28 ) 
for the squares. 

Two functions from [HI are 



or 



if x< -1, 
h 2 (x) = { sin(f sin 2 (f(l + x))) if - 1< x < 0, 

1 " if x > 0. 

if £ < -1, 

h^{x) = I cos (| sin 2 (fa;)) if - 1 < x < 0, 



if x > 0. 
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Figure 12. The figure shows several monotone increasing functions h(x) 
in [—1,0]. It is, in fact, showing left half of the window functions considered 
in this study. 

We also find another profile function following [H]. Here we consider a monotone in- 
creasing function (denote it as mallat on the legend of the figures) in [—1,0] 

{0 if x < -1 

sinff (l + sin(§ sin (f sin f^^))))) if - 1< x < 
1 X if x > 0. 

It is to note that if we consider two sine functions or four sine functions inside sin(7r/4(l + 
f(x)), then the decay in the frequency domain is not as fast as of considering three sine 
functions defined by h^x), so we decide to use h^x) as a representative for this class of 
profile functions. 

Here to find a good window function which has fast decay in frequency domain, we 
try several window functions g(x), where g(x) is defined from g as given above. It can 
be observed that the windows with profile function hi(x), and the window with profile 
function h±(x) are smoother than the other windows and have faster decay all over the 
frequency domain. 

We also define a dilated variant of the profile hi(x), effectively throwing out some of 
the zeros and ones on the outside of [0, 1]. We define the dilated profile function by 

h 5 (x) = /ii(l/2 + d(x — 1/2)), for some suitable < d < 1. 



In Figure 13, we compare several window functions to find out which one has faster 
decay in the frequency domain. Analyzing the graphs, we notice that the window func- 
tions defined from monotone increasing functions h 2 (x) and h^(x) have slow decay in the 
frequency domain. All other windows have quite similar decay in the frequency domain. 
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FIGURE 13. The figure shows several \g\ in a logarithmic scale. We 
consider the window function with profile function h\(x) with c = 
[1,1.3,1.5,1.7,1.9], h 5 (x) with c = 1,1.5, and d = 0.9, as well as con- 
sidering h 2 (x), h 3 (x) and h±(x). 

We observe that the window function formed by the monotone increasing function h±(x) 
and the stretched variant with d = 0.9 and c = 1.5 behave very well when £/7r is in the 
range of 20 to 60, a region that one might still expect to be relevant, with relative size of 
the Fourier coefficients between 10~ 4 to 10 -8 . 
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